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OO ' Abstract 

o 

' We address the exact resolution of a MINLP model where resources can be activated 

^SJ , in order to satisfy a demand (a partitioning constraint) while minimizing total cost. Cost 

functions are convex latency functions plus a fixed activation cost. A branch and bound algo- 
rithm is devised, featuring three important characteristics. First, the lower bound (therefore 



each subproblem) can be computed in 0(n log n). Second, to break symmetries resulting in 
, improved efficiency, the branching scheme is n-ary (instead of the "classical" binary). Third, 

a very affective heuristic is used to compute a good upper bound at the root node of the 
enumeration tree. All three features lead to a successful comparison against CPLEX MIPQ, 
which is the fastest among several commercial and open-source solvers: computational results 
showing this fact are provided. 
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1 Introduction 



— 1 ■ 

^ ■ Congestion phenomena arise in many real world problems. It is therefore a key issue in modeling 
O ! techniques to address those features. Latency functions (see [7]) are a useful tool to model 
nonlinear increase of the cost of a resource with respect to its usage: as an example think of 
traffic jams. Traffic problems are widely studied applying the concept of equilibrium (see [5]). 
Studies about behaviors of (possibly) congested systems where users act autonomously, i.e. 
selfishly, have addressed the problem of quantifying the loss of efficiency compared to a centrally 
optimized situation (see [8]). There's also been efforts in the direction of pricing edges of networks 
in order to influence users behaviors (see [1]). 

Our model addresses optimization problems where there is (i) a single commodity demand 
^ ■ that has to be satisfied, (ii) a set of potential resources, (iii) a fixed activation cost for each 
. resource and (iv) the congestion heavily influences the cost of a resource. We aim to a socially 
desirable solution, that is, a global optimum. 

We denote by Q the set of the q available resources. The cost ji{xi) of using resource i £ Q 
at level G [0, 1] is given by: 

, . _ J if = . . 

^^^^i) \ci + xji{xi) ifO<Xi<l 

where /j(-) — the latency function — satisfies the following: 

Assumption 1 /,(•) is semi-convex with /i(0) = and its derivative ^^^f^ > for ^ > 0. 
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Note that, for any i €z Q, the function Ci+Xif{xi) is not convex in the domain Xi > 0, although it 
is for Xi > 0. Unless differently specified, we assume troughout the paper that latency functions 
fi{-) satisfy Assumption [1] 
The model is the following: 



z* = min ^ ayi + Xifi{xi) (2) 
ieQ 

s.t. Xi < yi, for alH G Q (3) 

^x, = l (4) 

xeM\ (5) 
ye{0,ir (6) 

where Xi indicates the fraction of demand allocated to resource i (level of usage of resource i) 
and yi is a binary variable indicating wether resource i is active (xj > ^ T/j = 1) or not 
(xj = =^ = 0). The above formulation is a mixed integer nonlinear program and, in 
particular, all the constraints are linear with the y variables restricted to be integer. We will 
often refer to (HHS]) as P . 



2 Complexity issues 

In this section we prove that P is, in general, a NV-haxd problem. In fact even when /(•) is a 
homogeneous linear function P contains PARTITION, as we prove in Section 12.21 Nevertheless, 
as discussed respectively in Section 12. l.!) 12.1.21 and I2.1.3| the problem is easy when the cost 
functions are identical, or there are null fixed costs (i.e., Cj = for all i G Q) or when /(•) is 
constant. 

2.1 Easy cases 

2.1.1 Identical cost functions 

We consider the case when all cost functions ji{xi) are of the same type, that is, for all i £ Q, 
li{xi) = c+Xif{xi) = 7(xi). To the best of our knowledge, these results, though straightforward, 
are not present in the literature. In this case, the problem reduces to finding the optimal number 
k < q oi active resources. 

Proposition 1 // the resources are all of the same type with cost functions 7i(-) = 7(-)? for all 
i G Q, and an optimal solution for P consists of using k resources, then there is one solution 
where equal fraction of demand is allocated to the k resources. 

Proof. Based on the convexity of Xifi{xi), for any A;-uple of nonnegative numbers xi, . • • , x/j, 
with Yli=i Xi = 1, we have: 

k-l(l)<Y.^{xi). (7) 
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The next natural question we need to ask is "What is an optimal number of such resources 
(i.e., the best value for k)?". It is easy to observe that with zero setup costs (c = 0) the cheapest 
solution consists of using the maximum number (i.e., q) of resources. 

A solution that uses k + 1 resources costs no more than a solution with k resources if and 
only if the following is true: 

If c = 0, by definition of /(•), the last inequality is always valid. Therefore, it is cost-effective to 
use another resource if the additional setup cost does not exceed the gain in the variable costs. 
The cost of using k resources is 

If we define g{9) = 9f{9), which is convex by definition of /(•), we have that 

dk^ ~ dk^\ ^\k)) ~ \ky 

Hence, there must be 1 < A;* < g such that 

F{1) > F{2) >...> F{k*) and F{k*) < F{k* + !)<...< F{q). (11) 

Therefore, since a binary search can be used to efficiently find the k*, the following proposition 
holds. 

Proposition 2 When the q resources are all of the same type with the cost functions = 'y{-), 
for all i G Q, the problem is solvable in 0{Clog{q)) time, where C is the maximum computational 
effort for calculating 7(^). ■ 

2.1.2 Null fixed costs 

When the fixed costs are all null, then the objective function reduces to a sum of convex functions 
Xifi(xi), which is also convex. Since there is no need of binary variables to deal with fixed costs, 
we are left to the problem of minimizing a convex function over a polyhedron and this can 
be achieved via Karush-Kuhn- Tucker conditions. The optimal solution of this problem is also 
referred to as a Nash equilibrium (see [2]), because the partial derivatives of the cost functions 
are all set to the same value. 

2.1.3 Constant latency functions 

The objective becomes a linear function when the fi{xi) = fi are given constans. In this case 
an optimal solution is simply x*. = y*. = 1 where i* = argminjcj + fi, i G Q} and x* = y * = 
for all i i*. 



> 



(10) 



3 



2.2 NP-hardness of P 

From the discussion in Section [2.11 the minimal case for which the problem is possibly MV-hard 
is when there are q resources fi{x) = hiX. 

Proposition 3 Prohle'm\^^ is J\fV -hard even when the latency functions fi{ ) are linear. 
Proof. We use a reduction from the well known binary AAP-hard problem 

Partition: given q nonnegative integers wi,W2, ■ ■ ■ , Wq, with Yli=i — find, if 
it exists, a subset 5 C {1, . . . ,q} such that Yli<=s '^i — \- 

First, note that when we select a set of resources, the following convex subproblem may be 
easily solved by applying the Karush-Kuhn- Tucker (KKT) optimality conditions. Let 5 C Q be 
the set of selected resources. Then Problem [2H6] restricts to 

z{S) = min I /, + hixj ■.Y,^i = l-Xi£W\_\. (12) 
Vies ies J 

Restricted problem [12] is a convex optimization problem. Let us construct the following La- 
grangean function (for simplicity the constant term X^jgg fi has been omitted in the objective 
function below): 

L{x, //, A) = ^ bixf + A 1 - ^ Xi + ^ fiiXi 

The KKT conditions are (the symbol * indicates the optimal solution) 

V,L(x*,A*,/x*) = 0, 
^x, = 1 

i&S 

T 

fi* x* = 

X* > Og 



We then obtain 

A* 



and choosing x* ^ ^ fn = Vi G 5, we have 

. A* 1 



' 26. 26, Z,es i 

These values satisfy the KKT conditions and are therefore a global optimum for Problem 1121 

Consider now an instance of Problem [^HHl with fi = Wi and hi = ^ (x)' where the lOj's 
and W are those of the Partition instance. Considering the above Restricted formulation and 
given the set of active resources S, the optimal solution value can be rewritten as 

z{S) = min < u^j + > 
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Observe now that the function ^ + | has a minimum at ^* = -^/a, therefore z{S) is minimized 
when X^jg5 = y'^^^ ~ Hence Problem (HI)-® can be rewritten as 

min z(S) (13) 

5C{l,...,g} 

which is equivalent to finding 5 C {1, . . . ,q} such that X^jg5 Wi = This is, in turn, equivalent 
to Partition. This completes the proof. ■ 

3 Dual Bounds 

In this Section we will study the Lagrangean relaxation of the problem, its Continuous relaxation 
and the relation between them, in order to efficiently compute a dual bound. For the development 
to follow, we first address the case when the set S" C Q of selected resources (i.e. active facilities) 
is given. In this case, the fixed costs are known and variables Hi, i £ Q, vanish from (PHHI)- The 
resulting problem can be written as: 



ies lies ies } 



(14) 



and, since its objective function is convex, it may solved by applying Karush-Kuhn- Tucker 
optimality (KKT) conditions. 

Omitting the constant term ^jg^Q, the Lagrangean function 



L{x, ^, A) = X] Xifii^i) + 

i&s V ies J ies 



The KKT conditions, for the triple (x*,/x*, A*), are 

= /,(x*) + <//(<)- A* -M* = for allied 

E< = 1 

fi*x* = for all i G 5 
X* > for aWie S 
H* > for sllie S 

From X* > we have that fi* = for all i £ S. Then, 

x*=9-\y) (15) 

where gi{C) = /i(C) + C/i(C) — the derivative of C/i(C) — is a monotonically increasing functioifl 
and A* can be computed, independently of the x*'s, by 

Y^g^\x*) = l. (16) 



^Here we neglect the case constant fi. 
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Depending on the functions /i(-), the value of A* and x* can be computed in closed form or via 
numerical methods. For instance, if fi{C) = hC^, we have that: 



p I 

A* = I — and 



3.1 Lagrangean relaxation 

We now use Lagrangean relaxation to obtain a lower bound on z* , the optimal solution value of 
Problem (HH). 

Relaxing the activation constraints ^ using nonnegative Lagrangean multipliers Ki,i = 
1, . . . , IQI, we obtain the following problem: 

zlrp{k) = min J2ieQ i^i ~ '^i)yi + ^ifi^i) + f^i^i (17) 

s.t. ZieQXi = ^ (18) 

x€Rl (19) 

ye{0,l}i (20) 

Problem ()17H20p is a relaxation of the original Problem for any k >0q, and it is decomposable 
since optimal values for the y variables are independent of the values of the x. The optimal 
solution values for the y variables are the following: 

The remaining convex program, on the x variables, is: 

zlrp'{k) = min ^.gg Xif{xi) + KiXi (21) 

s-t. E.6Q^. = 1 (22) 
xeR% (23) 

and therefore zlrp{k) = zlrpi{k) + Y^^^q {ci — Ki)y*. Problem ([2TH23|) . is a convex program 
and can be optimally solved via KKT conditions. Using multiplier A G M for the covering 
constraint (j22p and multipliers fi £ for nonnegativity constraints (j23p , we obtain the following 
Lagrangean function: 

Lk(A, /i) = min ^ {xifi{xi) + KjXj - ^Xi) + A 1 - ^ Xj . 

ieQ \ ieQ J 

Thus KKT conditions are 

^h^^ll^ = g^(x*) + - X* - fi* = foralHEQ (24) 

OXi 

E< = 1 (25) 

H*x* = for alH G Q (26) 
X* > for all i G Q (27) 
IJ* > for alH G Q (28) 
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Now note that if i G 5 C Q then x* > and /i* = 0; since, from Assumption [H the first 
derivative is positive for poitive arguments, we have that gi{x*) = X* — Ki > 0, resulting in 

X* > Ki, for all i e S. (29) 

On the other hand, considering those i ^ S we have that x* = =^ /i| > =^ dii^i) = X*+fi* — Ki; 
the latter implication together with the fact that gi{0) = lead to n* = Ki — X* > 0. This implies 

X* < Ki, for ah i ^ S. (30) 

Equations [29] and [30] induce a particular structure for a solution to be optimal for ()21H23p , 
which is the following: 

Kl < K2 < ■ . . < Kh < X* < Kh+l < ■ . . < Kq. (31) 

' ^ ' ^ V ' 

S Q\S 
Furthermore from [25] and from the above discussion we obtain 

x* = gi\X*-K,) (32) 

and since we have X^jg^x* = 1 = X^ig5 5j"^(A* — Ki), the expression for A* is 

X* = G{Ki,...,Kh). (33) 

where G{-) can be computed in closed form (where possible) or via numerical methods. 

The above results suggest an efficient algorithm to compute the optimal solution of the 
relaxed problem (f2TH23]) : order the Ki in non-decreasing sense, compute A* using[33]such thatlST] 
is satisfied, compute the optimal variables for z G using [35] and set x* = for all i ^ S. We 
will refer to this algorithm as to the Ordering Algorithm. 

3.2 Continuous relaxation 

In the continuous relaxation of out original problem we relax the integrality constraint on the y 
variables. The problem becomes: 

w = m] 



^ ^ CiUi ~\~ Xifii^Xi) 


(34) 


i&Q 




Xi < Ui, for all i £ Q 


(35) 


= 1 


(36) 






X G 


(37) 


< < 1, for alii 


(38) 



It is easy to see that in the optimal solution of this problem x* = y* for all i G Q and thus we 
can drop constraints 1351 and 1381 and replace yi with Xj in the objective function. The problem is 
now the following: 

w = min ^ axi + Xifi{xi) (39) 

s.t. ^xi = l (40) 

X G Mi (41) 
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Note that setting q = Kj makes the Problem (j39H^T]) equal to Problem ([2TH23]) . that is 
w = zlrp'{c). By convexity (the proof is an extension of a well known result in LP, see [6] 
and [3]) we have that w = maxfc{2£,^p/ (A;)} and therefore we have the optimal multipliers for 
the optimal dual solution, that is k* = c. Another consequence of this observation is that the 
Continuous relaxation (whose optimal values equals the optimal lagrangean dual) is solvable 
using the Ordering Algorithm described in the previous Section. 

4 Primal Bound 

In this Section we describe a heuristic method that exploits the theory developed in the previous 
sections. The heuristic starts from the primal solution value induced by the optimal dual solution 
point (activated resources in the dual solution) and then performs a local search to improve it. 

Compute the set S* C Q of activated resources — and therefore A* — as in the Ordering Algo- 
rithm; then use Equation [15] to compute values of the activated variables {i G S) and set the 
others to zero [xi = 0,i ^ S). The current objective function value (primal) is: 



Now perform a local (greedy) search. Select — with a suitable criterion — an activated resource u 
to remove from S and compute z{S \ {j}): if this solution is better than the previous, update 
the solution (S := S \ {j})- Now select a non-activated v resource to include in S and compute 
z{S U {j}): if this solution is better than the previous, update the solution [S := S U {j}). 

The choice criterion for the local search can be tailored on the particular fi{xi). Anyhow a 
criterion which can be always applied is the one of choosing i £ S as argmaxjcj} and i ^ S as 
argminjcj}. As a example if fi{xi) = biX^, a criterion can be to compute the argmax of c,, of 
Ci + bi, and of Cj -|- '^^yfbl and to pick the argmax of them. 

5 Exact algorithm 

In this Section we devise an exact branch-and-bound algorithm for dSHl]) based on the develop- 
ments of the previous sections. 

5.1 Solution strategy 

In the enumeration tree, each node u represents a subproblem that is defined by {i) a set of 
active (ON) resources T, that is T = {z : = 1; i G {1, . . . ,g}} C Q, (ii) a set of inactive 
(OFF) resources F {F O Q, F n T = 0, F = {1 < i< q}) that must be OFF and {iii) the set 
Q \ {T U F} of resources that are not yet fixed to ON or OFF. If a resource is OFF, we consider 
unavailable for that specific subproblem. We indicate by Q(i/) = Q \ F the set of available 
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resources at node (subproblem) u. The generic subproblem may be formulated as follows: 
min ^ {ciyi + Xifi{xi)) + ^^{ci + Xifi{xi)) 

ieQ{u)\T i€T 

s.t. Xi < Ui for all i ^ Q{v)\T 

xi = i 

Xj > for all i € Qiv) 
Vi £ {0, 1} for all i £ Q{v) \ T 

Note that also for i £ T, the x'^s are still to be decided. Applying the same procedure described 
in Section 13.21 it is easy to see that the optimal dual value for a generic subproblem is given by 
the solution of: 

min ^ {axi + Xifi{xi)) + ^ (q + Xifi{xi)) (42) 

i€Q{u)\T i£T 

s.t. ^ Xi = 1 (43) 

i&Q{v) 

> for all i£Q{v) (44) 

This problem can be efficiently solved using the Ordering Algorithm: just set q = for all i £ T, 
because they must be in the solution since we pay their activation cost: the deeper is a node in 
the enumeration tree, the smaller is the problem to solve at that node, since the resources which 
are fixed OFF vanish from the problem. 

The upper bound at the root node is given by the heuristic method described in Section [H 

5.2 Branching strategy 

Two resources u £ Q and v £ Q are said identical copies iicu = and /^(x) = fvix), V x G [0, 1] 
(that resource is also said repeated). The number of copies of each resource (which is at least 
1), can be found by suitably preprocessing the instance. 

Branching here is to decide wether a resource i is used and in how many copies. Clearly, 
at node of the enumeration tree, we branch on resources i £ Q{i^) T. A possible branching 
rule — which turned out to be very effective — is to take i £ Q{i') \ T such that i = argmaxjci}. 

We adopt an n-ary branching scheme. If the variable to be branched corresponds to a resource 
i va. n copies, we generate n + 1 nodes (subproblems): in the l-th subproblem, / = 0, 1, . . . ,n, I 
copies of resource i are ON and n — I are off. Note that, with the n-ary branching, compared 
to a binary branching strategy (that is, not considering the impact of repeated resources on the 
solution), we only generate n + 1 subproblems instead of 2". 

6 Computational experiments 

The tests have been run on a PIV 3.2 GHz 1GB of RAM under Windows XP; our algorithms 
have been coded in C-I--I-. The computational experiments have been performed on the case 
where fi{xi) = biXi, for this case being the minimal NP-hard case and to be able to compare it 
with CPLEX MIPQ solver. 
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6.1 Design of instances 

We say that an instance is non- dominated if no resource is strictly preferable to any other, that 
is no i,j exists such that q + 6j < cj + bj and q < cj. To this purpose we impose the following 
relations: 

ci > C2 > • • • > Cg (45) 
61 < 62 < • • • < 6g (46) 
Ci + bi> ci for all i £ Q. (47) 

We considered two classes of instances. The first is referred to as base class and it has the 
following properties: (i) {bi ^ bj) and (cj 7^ Cj), for all 1 < i < j < q; (ii) = 6j + 1 for all 
1 < i < q and (in) bi = for all 1 < i < q. Instances in this class are deterministically 

generated and there are no repeated resources. In the second class (the random class), instances 
are randomly generated and repeated resources are allowed. We imposed a computation time 
limit of 2 hours. Base class and random class instances are denoted, respectively, with a "b" 
and a "r" in the instances' name. We considered up to 1300 binary variables for the base class 
and up to 100 binary variables for the random class (due to the time limit constraint): for each 
number of variables in this class we generated 20 random instances. 

6.2 Results and analysis 

Experiments results are summarized in Table [TJ Our algorithm outperforms CPLEX over all the 
instances classes and dimensions, except for very small problems (where the CPU time is very 
small indeed). 



instance 




Our algo 






CPLEX 10.0 






CPU time 


nodes # 


iter. # 


CPU time 


nodes # 


iter. # 


b200 


2.75 sec. 


2229 


1114 


1.97 


4608 


8112 


b400 


17.13 sec 


10897 


5448 


27.55 


42686 


69655 


b600 


66.89 sec 


34749 


17374 


184.95 sec. 


192012 


304124 


blOOO 


515.45 sec 


192591 


96295 


2492.58 sec. 


1486308 


2290007 


bllOO 


793 sec. 


274897 


137493 


4898.39 sec. 


2424233 


3700426 


bl300 


1774 sec. 


529275 


264636 


>7200 sec. 


> 4000000 


> 7000000 


r25 


0.120 sec. 


106.2 


36.2 


0.060 sec. 


352.8 


489.6 


r50 


0.450 sec. 


800.8 


134.6 


*106.370 sec. 


*731558.5 


*874274.3 


r75 


0.76 sec. 


1015.91 


131.4 


*227.27 sec. 


*1668124.9 


*1785904.7 


rlOO 


5.3 sec. 


7774.45 


633.1 


*821.35 sec. 


*4575207.7 


*5277239.6 



Table 1: Random instances (r) and base class instances (b); an asterisk means that CPLEX was stopped 
due to memory or CPU time limit (2 hrs). Data for random instances are average. 

We first highlight the excellent performance of the heuristic procedure: in all the instances 
it found the optimal solution in a considerably small CPU time (always within seconds). As 
it can be seen from the dimension of the instances we tested, the more difficult ones are the 
random instances; in many cases CPLEX was stopped because of time limit or because of 
memory usage limit. This never happened for our algorithm. (The average values in those cases 
is computed considering the limit time or the time at that CPLEX stopped.) The reason why 
random instances are more critical is that repeated resources are allowed. Our multidimensional 
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Figure 1: CPU time comparison for random instances. Om' algorithm outperforms CPLEX. 

branching strategy turned out to be very effective, while CPLEX does not have any mechanism 
to avoid considering nodes which are exactly the replication of other nodes (when branching 
on repeated nodes). Figure [T] highlights the performance of our algorithm versus CPLEX: the 
bigger the instances the better is our algorithm. 

7 Conclusions 

In this paper we studied the problem of allocating a demand when the resources' costs are 
modeled with general convex functions. We devised a general solution approach based on a 
branch-and-bound scheme. The upper bound at the root node is given by a heuristic method 
whose perfomance is excellent (it found the optimum in all the tested instances). The lower bound 
at the leafs is a dual bound that is computed solving the continuous relaxation of the problem 
(their optimal values are identical) and for that task we devised an effective exact algorithm (the 
Ordering Algorithm). Computational experiments show that we outperform CPLEX and that 
our particular multidimensional branching strategy is very incisive. 
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